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We apply algorithms based on Lieb-Robinson bounds to simulate time-dependent and thermal 
quantities in quantum systems. For time-dependent systems, we modify a previous mapping to 
quantum circuits to significantly reduce the computer resources required. This modification is based 
on a principle of "observing" the system outside the light-cone. We apply this method to study 
spin relaxation in systems started out of equilibrium with initial conditions that give rise to very 
rapid entanglement growth. We also show that it is possible to approximate time evolution under 
a local Hamiltonian by a quantum circuit whose light-cone naturally matches the Lieb-Robinson 
velocity. Asymptotically, these modified methods allow a doubling of the system size that one can 
obtain compared to direct simulation. We then consider a different problem of thermal properties of 
disordered spin chains and use quantum belief propagation to average over different configurations. 
We test this algorithm on one dimensional systems with mixed ferromagnetic and anti-ferromagnetic 
bonds, where we can compare to quantum Monte Carlo, and then we apply it to the study of 
disordered, frustrated spin systems. 

>h : I. INTRODUCTION 

Matrix product and density-matrix renormalization group methods provide one of the most powerful ways of 
Ch ' simulating one-dimensional quantum systems. In addition to ground state properties they have been extended 
to thermal and open systems [2,1 and dynamical problems 3]. The reason for the success of these algorithms is that 
in many cases the appropriate quantum state can be well-approximated by a matrix product state, giving a very 
compact representation of the state of the system. 

In some cases, we even have theorems that quantify the accuracy of matrix product states. For ground states of 
local quantum systems with a gap, the ability to represent ground states as matrix product states follows from bounds 
on the entanglement entropies [J| . Related results are available for thermal systems [a, @] and for non-equilibrium states 
■ obtained by starting with a factorized state and evolving under a local Hamiltonian for a time t; in the first case 
^-H the bond dimension needed to get a good matrix product approximation to the desired state scales exponentially 
fN| ' in /3, while in the second case it scales exponentially in i[7.]. The two results @, 0| are constructive proofs, which 
give an algorithm to find the matrix product state. All of these constructive proofs rely heavily on Lieb-Robinson 
Q ' bounds[i,[i,[ll[il|. 

_ In contrast to these constructive proofs, the matrix product algorithms used in practice are variational: they involve 
. optimizing over different matrix product states to find the best one. This works well because in many practical cases 
^ ■ the entanglement grows much more slowly than the upper bounds set by theory. For example, systems described by 
• i-H , conformal field theory have an entanglement entropy growing only logarithmically with system sizep^. while the area 
' law bound[4| gives no useful result for these systems due to the absence of a gap. Similarly, for many initial conditions, 
^ ' evolution under a local Hamiltonian gives an entanglement entropy growing only logarithmically in timep^ [l^ . [Tsj 
- - ' while the theoretical upper bound gives an entanglement entropy growing linearly in time [1^1 ■ 

In this paper, we argue that there are many situations in which algorithms based on Lieb-Robinson bounds are 
the best technique. We first look at the case of time evolution in systems out of equilibrium. Here, there are initial 
conditions for which the entanglement entropy is known to grow linearly in time[l7|, in accordance with conformal 
field theory predictions fl8| . Roughly speaking, logarithmic entropy growth tends to occur in cases where we can 
divide the chain into a small number of subchains such that the initial state is an eigenstate of the Hamiltonian on 
each subchain; for example, starting an XXZ spin chain in a state in which all the spins on the left half of the chain 
are up and all the spins on the right half of the chain are down leads to a logarithmic entropy growthfisj. On the 
other hand, the linear entropy growth tends to occur in cases where the initial state differs from an eigenstate of 
the Hamiltonian on every subsystem of the full spin chain. For example, starting an XXZ spin chain in an initial 
condition in which the spins alternate between up and down (a Neel state, the ground state when the Ising term is 
the only term in the Hamiltonian) leads to a linear entropy growth [l7j. 

For system with linear entropy growth, matrix product methods will require a bond dimension growing exponentially 
in time to obtain accurate results. At this point, both variational matrix product and constructive, Lieb-Robinson- 
based methods require resources growing exponentially in time. The question, then, is how to obtain the smallest 
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exponential. To some extent, the Lieb-Robinson based methods (such as [tJ and the methods below) are "worst-case": 
the bond dimension depends on theoretical upper bounds for arbitrary local Hamiltonians, while the matrix product 
methods can adaptively find better representations. On the other hand, there are some disadvantages to matrix 
product methods. To get a rough idea of the resources required, let us consider a system of TV spins, each of spin-1/2, 
with a local Hamiltonian. The simplest algorithm to simulate this system for a time t involves writing the initial 
condition down in some basis and then directly simulating it (we discuss below different methods for doing this), 
requiring resources scaling as 2^t using sparse matrix methods. A matrix product algorithm can avoid truncation 
error for this system by using matrices with bond dimension 2^/^19, 20]. However, the algorithm must then perform 
singular value decompositions and eigenvalue calculations, taking a time which scales as the cube of these matrices, 
and hence of order 2^^/^, which is slower. Of course, the matrix product methods are really only useful if they are 
able to represent the system with a smaller bond dimension. In this case, though, even if such a representation exists, 
the algorithm must find it, and this can pose a problem. The algorithm for simulation of non-equilibrium systems 
depends on breaking the time evolution into a series of small Trotter steps; if the Trotter steps are too long this can 
lead to Trotter error, while if they are too short, the truncation error can grow rapidly: even if there is a good state 
the algorithm may not find itTB*!. 

We can use this scaling of the difficulty with N to get an idea of the scaling of computation effort with time, using 
a Lieb-Robinson bound on the group velocity, vlr- Suppose we wish to compute the expectation value of a local 
observable, such as a spin on a site, at a time t/, starting from a factorized state at time t = 0. The number of spins 
in the past light-cone of this spin is 2vLjit, and the Lieb-Robinson bounds imply that the effect of spins outside the 
light-cone is exponetially small. Thus, it suffices to simulate only the dynamics of the N = 2vLRtf + 0(log(e)) spins 
closest to the given spin in order to compute the expectation value to an accuracy e. This requires an effort scaling 
exponentially in time as . Similarly, if the entanglement entropy grows linearly in time, the matrix product 

methods also require an effort scaling exponentially in time. 

Our main result in this paper is the light-cone quantum circuit algorithm, an application of the Lieb-Robinson 
bound that allows us to simulate the evolution of local observables with resources growing asymptotically as only 
Nt2'"^'^, by some statistical sampling. This allows twice as large systems as the direct method. We analyze the entropy 
growth in these systems and argue that matrix product methods are also less efficient for long time simulation. We 
apply the light-cone quantum circuit algorithm then to the problem of spin relaxation in spin chains started in the 
Neel state. The physical idea behind the light-cone quantum circuit algorithm is as follows: to find the state of a 
given spin at a time t / , we only have to track the dynamics within the past light-cone of the spin. For times t close 
to zero, the past light-cone includes roughly 2vLRtf spins, but at these early times the entanglement is small and 
hence the computational effort should be less. For times t close to t/, the past light-cone includes few spins and hence 
should be easier to simulate. 

The paper is organized as follows. We first derive the light-cone quantum circuit algorithm. We then apply it to 
spin relaxation, and study oscillations of the central spin, decay of the envelope of the oscillations, and also seemingly 
random oscillations of the central spin in chains where boundary effects become important. We then derive a related 
quantum circuit method, the corner transfer quantum circuit which may be useful for studying the evolution of global 
observables in highly entangled non-equilibrium states. 

We then turn to a different problem, presenting one other application of Lieb-Robinson methods, using the quantum 
belief propagation algorithm[22'| to study thermal states in disordered systems. The quantum belief propagation 
algorithm explicitly constructs a matrix product state for a thermal quantum system. While it manipulates operators, 
rather than states, and thus can be computationally expensive, it has other advantages. It has no Trotter error, making 
it fast and accurate at high temperatures: it can obtain quantities such as the susceptibility peak to higher accuracy 
using fewer resources than methods such as transfer matrix DMRG 23], although at low temperatures it breaks down, 
with the resources required scaling exponentially with the temperature. In this case, the exponential scaling with 
the temperature is again related to a linear relationship between a time scale, in this case /3 = l/T, and a length 
scale. It can be applied to random systems, where transfer matrix DMRG cannot be used because of a lack of 
translation invariance. A good test of variational matrix product [2] methods on this kind of system is lacking, so 
we cannot compare here. We apply the quantum belief propagation algorithm to two random systems, one without 
frustration where we can compare to quantum Monte Carlo and one with frustration where Monte Carlo methods are 
not applicable. 

II. QUANTUM CIRCUIT METHODS 

In this section we present the various quantum circuit methods. We begin by reviewing previous work, and then 
derive the light-cone quantum circuit algorithm, and apply it to a problem of spin relaxation. We then use previous 
results on the entanglement entropy growth to estimate the computational resources required for different approaches 
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to this problem, and finally we present the corner transfer quantum circuit method, an extension which allows access 
to global quantities. 



A. Background 

To understand our algorithm, we first review the ideas in fj^ which gives a construction of a matrix product 
operator approximation to the time evolution operator, exp(— using resources exponential in t. We consider a 
local Hamiltonian 

H = J2h, (1) 

i 

where each i acts on sites i,i + 1. This Hamiltonian obeys a Lieb- Robinson bound: given any operator O which has 
support on set of sites X, the operator exp{iHt)0 exp{—iHt) can be written, with exponentially small error, as an 
operator acting on the set of sites i within distance VLRt of X, where vlr is the Lieb- Robinson group velocity. 

To simulate the system for a time t, we divide the system into blocks of length Z, where I is slightly larger than 
2vLRt (the error in the approximation will be exponentially small in / — 2vLRt). We then let H = Hq + H' where Hq 
is the sum of the Hamiltonians on each block and H' is the Hamiltonian connecting the blocks: 

i<{k+l)l-l 

^o = E E (2) 

k i=kl + l 



H' = Y,hku (3) 



k 



where the sum ranges over integers k. 
We then write 

exp{-iHt) = (t exp[-i^ cxp{-iHot')H' cxp{iHat')]^ exp{~iHot), (4) 

where T denotes that the exponential is time-ordered. The operator exp(— ii/o^) is equal to the product Yik Uk, where 
Uk is a unitary operator acting on sites kl + 1, kl + 2, kl: 

i<(fe+l)i-l 

C/fc = exp(-i ^ hit). (5) 

i=kl + l 

Using the Lieb-Robinson bounds, we can approximate exp{^iHot')hki expliHot') by an operator hki{t'Y°'^ which 
has support on sites kl — 1/2 + l,...,kl + 1/2 — 1 for \t'\ < t and for I > VLRt- Then the operator 
Texp[— i Jq eicp{—iHQt')H' exp{iHQt')] can then be approximated by a product Yik where 

VI. =Texph^ / hu{tT% (6) 
Jo 

The key in this construction is that the intervals kl — 1/2+1, ...,kl + 1/2— 1 do not overlap for different k. This 
construction expresses the time evolution as a quantum circuit: 

expi-iHt)^Y[VkYlUk. (7) 

fc k 

The support of the operators Uk, Vk and the circuit is shown in Fig. [1] 

Precise error bounds can be given using Lieb-Robinson bounds for the error in Eq. ([7]). To get an error of order e 
in the propagator ([7]), we only need to take I = 2vLRt + 0{\og{N/e)). In what follows, we will not make detailed error 
estimates, since Lieb-Robinson error estimates are fairly simple and are by now standard in the literature; when we say 
that it suffices to take a length scale "of order" VLRt to obtain an approximation to a given local quantity, we mean that 
by taking the length scale VLRt + 0(log(iV/e)) we can obtain an error of order e in the state; when we are computing 
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FIG. 1: Support of the operators Uk,Vk and the quantum circuit for A'' = 12, I = 4. The smaU numbers at the bottom 
of operators Uk,Vk represent the bond dimensions required to represent these operators as matrix product operators; the 
maximum product of these across any bond is 16 = 4'^^ . 



expectation values of local quantities, to obtain an error of order e we need a length scale v^nt + 0(log(wL/jt/e)), so 
that the error bound does not depend on N in this case. 

Suppose we want to apply the quantum circuit procedure to compute the time evolution of some state ^'o. For 
simplicity, let ^0 be a factorized state (later we discuss the case where ^0 is a matrix product state both in this 
procedure and using our algorithm, and we find that using the idea of "observation" discussed below the case of 
matrix product state initial conditions presents no additional difficulty). The operator Uk is an operator acting on 21 
sites. Any such operator can be written as a matrix product operator with a bond dimension equal to 4^'/^ = 2^' jl^. 
This maximum bond dimension is achieved halfway across the interval of length 2/, while at any point a distance d 
from one of the ends of the interval the bond dimension is only 2"^, as shown in Fig. [1] 

The same holds for 14, and so the maximum product of bond dimensions across any bond is 2^' (a slightly worse 
estimate of 2"*' was found for this construction in since the fact that the bond dimension may vary with position 
was not taken into account). 

There are several problems, however, with implementing the above method in practice, and it is these problems 
which we overcome with our methods, the light-cone quantum circuit method and the corner transfer quantum circuit 
method. The first method is most appropriate for computing local quantities (such as a spin or energy expectation 
value) while the second method is most appropriate for finding a good global approximation to the ground state. 

The first problem is that operator equations of motion are computationally expensive in practice. To this end, we 
will modify the procedure to deal only with state vectors, rather than operators. The second problem is that the 
velocity of the quantum circuit does not obviously match the Lieb-Robinson velocity, in the following sense: given 
arbitrary operators Uk,Vk, supported as described above, the product life ^fc life ^fe ^an propagate information by 
a distance of 21 in each time step. Since I is roughly VLRt, this means that such a quantum circuit could have an 
effective velocity roughly 2vlb.- Of course, the operators Uk, Vk are not arbitrary operators, but still we would like to 
fix this problem; we will show how to do this with the corner transfer quantum circuit method below which also leads 
to improved estimates on the maximum matrix product state dimension needed. 

However, the real problem with this method is that it doesn't lead to any improvement over naive simulation when 
it comes to computing local observables. The main problem we consider in this section is the following: we start 
a spin chain at time t = in a factorized state ^^o, and then evolve under a local Hamiltonian for to a final time 
i/, at which point we wish to compute some local observable, such as S'f , the z-expectation value of spin i. By the 
Lieb-Robinson bounds, we can approximate Sf{t) by considering only a subchain of the full chain: wc consider only 
sites i — + I, where I is slightly larger than VLRtf. We then define V?' to be the appropriate factorized state 

on this subchain and evolve for a time tf using the Hamiltonian H' acting on the subchain. We then compute 
exp{iH'tf)Sf exp{—iH'tf)\'i''). Using sparse matrix methods to compute the time evolution of this requires 
a computational effort of order 12^'- , which scales as 2^"^''*^'. The quantum circuit method discussed above would also 
require simulations on intervals of length 21, and hence leads to no improvement when computing this local quantity. 
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FIG. 2: Support of the operators Hl, Hr, Hm,H'i^, H'^. 



B. Light-Cone Quantum Circuit Algorithm 

We now show how to reduce the computational effort to an amount of order 2'"^"-*f, allowing the time scale to 
be twice as large, by combining Lieb-Robinson bounds with statistical sampling. While we focus on this section on 
starting in a factorized state, we later discuss the case of starting in a matrix product state. 

We now derive our algorithm, which we call the light-cone quantum circuit method as it avoids keeping track of 
certain degrees of freedom outside the light-cone by making certain "observations" to reduce the computational effort. 
We first define a subchain of length 21 + 1 and an initial state on that subchain as above. We label the sites in the 
subchain by 0, We let H' be the Hamiltonian on the subchain and we write 

H' ^ Hl+Hr + Hb. (8) 

where TJ^ acts on the left half of the chain (sites —1), Hr acts on the right half of the chain, and the boundary 

Hamiltonian Hb acts on sites —1, 0, 1. Thus, [iJ/,, Hfj] — 0. We define Hm to act on the middle half of the chain: it 
is supported on sites —1/2, +1/2 (we pick / even for simplicity). See Fig. O 

Using the Lieb-Robinson bounds in the same way as in the quantum circuit method above we can approximate the 
time evolution for a time U = tf/2 by: 

'^'{tf/2) = exp(-ii7't//2)*' w exp(-ii?Afi//2) exp[-i{HM - HB)tj/2] exp[iHLt//2] exp[ii/iit//2]*'. (9) 

Note that the decomposition ([9]) is only good for times of order tf/2. We wish to compute the expectation value of 
Sq, or any other observable on site 0, at time tf. Again using the Lieb-Robinson bounds, the expectation value of 
this is approximately equal to 

(*'(ij/2)|exp(iFMi//2)Ocxp(-iJJMi//2)|*'(t//2)). (10) 

Combining Eqs. (|9TT0l) and our use of the Lieb-Robinson bounds to approximate the expectation value of O on the 
full chain by its expectation value on the subchain, we have 

(^'o|0(t)|*o) « (*|0|*), (11) 

where 

* ^ exp(-ii?A/i//2)4'' w exp(-ii7Mi/)exp[-i(ffM - ffi3)t//2] exp[ii/Lt//2] exp[ii/i?t//2]*'. (12) 

The operator Hm — Hb is a sum of two operators, H'j^ and ff^, where H'^ acts on sites —1/2, —1 and iJ^ acts on 
sites +1, .... +1/2. Therefore, 

^ = exp(-ii/Afi//2)*' « ex.p{-iHMtf)(cxp[iH'^tf/2] exp[-iHLtf /2]j (ex.p[iH'Btf /2] exp[-^i^i^^//2])^''. (13) 

Eq. (fT3|) is not yet useful computationally, since it will require an effort of order ltf2^^ to compute the evolution of 
the state ^' . We now describe the light-cone quantum circuit method to compute the expectation value: first, write 



(14) 
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where "^l, are states on the left and right half of the chain, and is a state on the center site. We compute the 
states 

■^'l = (exp[^i^i^//2]exp[-^i^L^//2])*L, (15) 

and 

^'j, = {exp[iH'jitf/2] exp[~iHRtf/2]^^R, (16) 

which requires an effort of order ltf2K Next, we introduce a complete orthonormal basis of states on the sites 
— /, ~l/2 — 1, which we label 0l(Q!), and another a complete basis of states on sites 1/2 + 1, +/ labelled (j)R{a). 
Then we decompose ^E*'^ and ^'^ as: 

vI'i=^A(a)0L(a)§5a(a), (17) 

a 

where CL(ck) is some normalized state on sites — Z/2,...,— 1 and ^^(a) is some normalized state on sites +l,...,+^/2. 
The states (/>l(q:) need not be eigenvectors of any reduced density matrix and the states Cl{<x) need not be orthogonal 
to each other, as Eq. (|17p is not a Schmidt decomposition. Thus, from Eqs. (|13ll5ll6ll7p . 

(*|0|*) =J2T. \A{aL)\'\A{aR)\'E{0, aL,aR), (18) 

where 

E{0, aL,aR) = {^L{aL) «) *c » ^R{aR)\exp{iHMtf)0 exp{~iHMtf)\^L{aL) «) *c ® Cfl(ai?))- (19) 

Eq. is at the heart of the light-cone quantum circuit approach. Numerically we proceed as follows: first, we 
compute |A(ai)p and |A(afl)p for all aL,aR. This requires an effort of order 2'. We then do a statistical sampling: 
we randomly pick an aL and an aR according to the probability distributions |A(ai,)p and |yl(ai?)p, and compute 
the average E{0, a^, a/?). We repeat this procedure many times, to average over different choices of a^, aR. 

The computational effort required then scales as only 12\ or roughly ^2"^^* as claimed. Asymptotically this allows 
double the time. The time scales linearly in the number of iterations of statistical sampling, which we denote Na- 
However, if O has bounded operator norm, then E{0, a^, aR) has bounded moments, and so by the central limit 
theorem, the number of iterations required still scales only polynomially in the error. 



C. Results on Non-Equilibrium Dynamics 

We now discuss results from this method, as well as some implementation details. We consider evolution under the 
XXZ Hamiltonian 

H = Y.(S!Sf+, + SfSf^, + A5f 5f+i). (20) 

i 

For A = 0, this problem can be mapped to free fermions by a Jordan- Wigner transformation and solved exactly. We 
use this as a check on our results later. 

We used as a starting point the Neel state, with spins alternating up and down, and we computed the time 
dependence of S'^(t) for the central spin. The main numerical effort is to compute the evolution of a state under 
a Hamiltonianm, which we did using a combination of short steps with a series method. For example, to compute 
exp(— i7?it)\['i we divide the time t into shorter intervals of time to, and compute exp(— ii/ito)^'L = "^l —itoH^'i'L — 
(iQi?£/2!)5'i -I- keeping a fixed number of terms in this series. We then repeat this procedure (t/to) times. To 
obtain negligible error for a chain of 20 sites with to = 1 required going to roughly 40-th order for |A| < 1, while for 
A = 2 slightly longer series were required. A more sophisticated way of doing the time evolution would be to build 
a tridiagonal Hamiltonian in the Krylov space spanned by ^l, H^'^ H^'^ ^ for some k, and then evolve exactly 
with this Hamiltonian[25]. 

Another important point of numerical simulation is the use of symmetries. We can choose the states (f'Lict) and 
4>R{a) to be eigenstates of total S^. Then, since the state ^P' is an eigenstate of total S^, the state ^L{ctL)'^'^c'^^R{ctR) 
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is also an eigenstate of total S" , which allows us to use symmetries when computing the evolution of the state. Since 
most of the numerical time is consumed statistically sampling E{0, a^, an), we build the sparse matrix for the 
Hamiltonian Hm in each spin sector once, before doing the sampling, and then run the sampling. 

It is also possible, although we did not implement it, to take into account reflection symmetry. Since and an 
are chosen independently, the state ^lI^l) ^ £,r{o.r) does not have reflection symmetry. However, if both the 

Hamiltonian Hm and the operator O — Sf have reflection symmetry, then it is useful to write 

CL(aL)®*c«)Cfl(ai?) ^ s{aL,aR) + ^ A{aL,aR), (21) 

where '^Si^A are symmetric and anti-symmetric states. Then, 

E{0,aL,aR)=^ {^s{aL,aR\exp{iHMt)Oexp{-iHMt)\)'i'siaL,aR)) (22) 
+ {^AiaL,aR\ exp{iHMt)0 eyip{-iHMt)\)^ A{aL, aR)) , 

and so we can statistically sample one of the two terms on the right-hand side of Eq. (p2)) . Note that on each iteration 
we randomly choose an and an aR and then randomly choose a term in Eq. (|22p. rather than repeatedly sampling 
Eq. 

As the algorithm is described above, the initial computation of the states ^'^(a) and ^'_R(a) depends on the final 
time. For each final time t, we have to compute a new set of states ^l{o() and ^R{a) and then do the statistical 
sampling. However, in fact, we can speed the algorithm up at a slight cost in accuracy: we fix a given tf and on each 
statistical sample we compute the state 

exp{-iHMtf)\^LiaL) (E) *c (23) 

to evaluate the expectation value in Eq. (fTg]) . We then act on this state with the operator ex-p{iHMSt) for some small 
St. We then use this new state to compute an approximation to the expectation value at time tf — St. We then act 
on that state with exp(—iHMSt) to compute an approximation to the expectation value at time tf — 2St, and so on. 
Since the computational cost of performing the time evolution of a state under a Hamiltonian is proportional to the 
time evolved, these additional steps are relatively cheap, for St << tf. There is a small cost in accuracy: in general, 
to compute expectation values at a time t, we can do the initial evolution for a time ti, and then evolve further for 
time t — ti. To make the effect of boundary conditions as small as possible, we would like to have both ti and t — ti 
as small as possible, which is why above we chose to evolve for a time ti — t//2. However, if St << tf, then we are 
not far away from the ideal choice of ti by initially evolving for time tf/2 and then evolving for time tf/2 — St. We 
followed this procedure in the numerical work below, with St — 0.25 and taking to be spaced with integer steps 
using to ^ 1. This accounts for some of the slight kinks in the curves after every integer value of t. 
The spin-wave velocity of Hamiltonian ((20|) for A < 1 is given by 

v,^ = {7r/2)sm{0)/9, (24) 

where cos(6') = A[23|- On the other hand, in the development above, we used a Lieb-Robinson velocity vlr, where 

vlr>Vsw (25) 

Using the Lieb-Robinson bounds, we showed that we could accurately simulate for a time t using length scales I = v^Rt. 
As mentioned above, we do not give precise error estimates, but it is not hard to give rigorous estimates of the error. 
However, the Lieb-Robinson bound is actually fairly conservative because of Eq. ()25p . and thus in practice length 
scales Vswt suffice to get good results. 

We begin by illustrating results for the XY chain, with A = 0. In Fig.[3l we illustrate results from exact simulations 
of the XY chain for various sizes. We consider chains with open boundary conditions with N = 35,51,101 and a 
chain with periodic boundary conditions with N — 36. For the chains with open boundary conditions, we plot the 
average of the central spin as a function of time, while for the periodic boundary conditions we plot the average of 
an arbitrarily chosen spin as a function of time. The large exact simulations are possible because this chain can be 
mapped to free fermions by a Jordan- Wigner transformation. Later we will present comparison of these results to 
light-cone quantum circuit results. For now, we discuss aspects of Fig. [3] which show the influence of boundary effects. 
In the range of times, the chain with N = 101 shows no effect of the boundary. There are oscillations of the spin with 
frequency uj roughly 2, with an envelope decaying as l/\/t so as 



(S^) ^ — cos(t^t-|-0o), 



(26) 
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FIG. 3: Time dependence of {S^} for the central spin as a function of time. Curves are A'' = 101 (black), A'^ = 51 (red), A'^ = 35 
(green) , and A'' = 36 (blue) . 



for Oq « (3/4)7r. Simulations on longer chains show that as long as A'' is less than t/2, the decaying oscillations of 
Eq. (p6)) continue to hold. In regard to the decay of oscillations found here numerically, it is interesting to note 

a similar power-law decay found for a different system of free bosons, where fluctuations about a maximal entropy 
state were proven to decay at least as fast as 1 /t^/^[2ll|. 

For N larger than t/2, the boundary conditions become important, as one can see in the curves for A'^ = 35 and 
N — 51, which start to show deviations from the A^ = 101 curve for times roughly 16 and 22 respectively. This is no 
surprise, since the distance between the central spin and the boundary is N/2, and Vsw = 1 for A = 0. We will see 
later for A =/= that in general boundary effects become important for A'^/2 = tvsw Interestingly, periodic boundary 
conditions offer no improvement, since the A^ = 36 periodic curve deviates at the same time as the N = 35 open 
curve. 

Another interesting effect is that once the boundary conditions become important, the expectation value shows 
wild oscillations, which no longer decrease in magnitude. This may be a consequence of the fact that the chain is 
integrable. It would be interesting to see for a non-integrable system whether such oscillations occur or not; the 
simplest statistical assumption for a non-integrable system is that the state at long times would be a random pure 
state satisfying the conservation laws of total and total energy, and thus the expectation value of {S^) for any i 
would show only exponentially small fluctuations about the average spin. 

We now consider the application of the light-cone quantum circuit method to this chain. We considered ^ = 18 and 
20, and did A^^f = 1000 iterations of statistical sampling, as shown in Fig. [H The results for exact simulations with 
A'' = 35 and A^ ~ 101 are also shown for comparison. We see that the light-cone quantum circuit method with I — 18 
is accurate over the same range of times as the exact simulation with A^ ~ 35, while the light-cone quantum circuit 
method with I — 20 improves on this result. By increasing I from 18 to 20, we increase the range of times by roughly 
2, while increasing A^ by 2 would increase the range of times by only 1. 

There are statistical fluctuations in the light-cone quantum circuit results in Fig. |31 due to random fluctuations in 
E(0,aL,Oiii) for different choices of aL,OiR- In Fig. [5] we plot the rms fluctuation in E{0,aL,aii) as a function of 
time, sampling this expectation value with the probability distribution |^(aL)p|yl(afl)p. The results in Fig.[4]are an 
average over Nu = 1000 samples, so the spread on each data point in Fig. H] is equal to 1/VTOOO times the fluctuation 
shown in Fig.[5l or roughly 0.005 in the worst case. A very interesting point is that for t < 5, the rms fluctuation is 
negligible, while for times t > 9, the rms fluctuations are up to their full value, with a sharp bend in the curve (plotted 
on a log scale) around < = 8 or i = 9. This is a consequence of the maximum spin-wave velocity: the influence of the 
regions on sites —I, —1/2 and +//2, +1 takes a time of order l/2vsw to reach the central region. There are sill 
some fluctuations for t < l/2vsw = 9, but they decay rapidly until they are negligible at very short time. 
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FIG. 4: Time dependence of (S^) for the central spin as a function of time for A = 0. Exact curves are A'' = 101 (black) and 
N = 35 (green). Light-cone quantum circuit curves are I = 18 (red) and I = 20 (blue). 




FIG. 5: RMS fluctuations in E{0, ul, an) as described in text as a function of time for I = 18. 



We then applied the light-cone quantum circuit method to chains with A — 0.5, 1, 2 as shown in Figs. I6I7I8I For 
comparison, we show an exact simulation of a chain with N = 20. For larger A = 0.5, 1 we still see decaying 
oscillations, but the decay is much more rapid than for A = 0. The decay of the envelope is, very roughly, t~^-'^^ 
for A = 0.5. For A = 2, no oscillations are seen. The effects of statistical noise are much more noticeable, since 
the magnitude of the spin is much less. All simulations were done with Nu = 1000 statistical samples, except the 
simulation with A = 0.5, 1 = 22 has only Na = 250 samples, and the simulations with A = 1, ^ = 16 and A = 1, ^ = 18 
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FIG. 6: Time dependence of (S^) for the central spin as a function of time for A — 0.5. Exact curve is = 20 (black). 
Light-cone quantum circuit curves are / = 18 (red), I = 20 (blue), and I = 22 (green). 



had Nit = 10000 and Nu — 3000 samples respectively. The spin-wave velocity is larger for these chains, so the 
simulation breaks down at an earher time than for A = 0; we again see that the simulations work for 

t < N/{2v,^), (27) 

or 

t < l/Vsn,, (28) 

in the exact and Hght-cone methods respectively. 



D. Entanglement Entropy 

From the predictions in [isj and the numerical work in [l3], we can compare the difficulty of doing a similar 
calculation doing time-dependent DMRG or time evolving bond decimation. The entanglement entropy is predicted 
to grow linearly in time. In [l3], the pref actor of the numerical growth was determined in a quench of an XXZ chain 
from an Ising coupling with A = Aq to an Ising coupling with strength Ai. The prefactor depended on Aq and was 
largest in the case of Aq = oo, the case we have considered here where the initial condition is a Neel state. There, 
in a quench to Aq = 0, the entanglement entropy (with logs taken to base 2) was observed to grow a little faster 
than O.Qvswt- This implies that the minimum size of the bond dimension needed in a matrix product state is at least 
2°-^'"=™*, which requires an effort going as 2^-*"=™*. In contrast, the numerical effort required to directly simulate a 
chain of length TV scales as 2^. To get accurate results for a time of order t, requires N = 2vswt, or an effort 2^"""*. 
Using the techniques in the present paper, the effort can be reduced to of order 2'"="'*. 

When comparing to IitI . note that there is a difference in normalization of the Hamiltonian by a a factor of 4, 
but since the results in [l7| are expressed in terms of the spin-wave velocity multiplied by time, the time axis is the 
same for A = 0. For A > 0, the time axis in [l3] differs from the time axis here, since [l^ multiplies the time by the 
relevant spin- wave velocity which is greater than unity. For A = 1, the spin- wave velocity using our normalization is 
equal to tt/2 > 1, and for A = 0.5 the spin-wave velocity is equal to 3-\/3/4, so the times in the present paper should 
be multiplied by such a factor greater than unity when comparing to [I3 . 

In [l3|, other initial conditions were considered, other than just the Neel state. These initial conditions were chosen 
to be the ground state of various other XXZ Hamiltonians. As the initial A was reduced, the entropy growth was 



11 



0.5 




-0.1 







5 



10 



15 



Time 



FIG. 7: Time dependence of (S'^) for the central spin as a function of time for A = 1. Exact curve is A'^ = 20 (black). Light-cone 
quantum circuit curves are I — 16 (red), / — 18 (green), and I = 20 (blue). 



found to still be linear in time, but with a smaller prefactor. Such simulations could still be carried out with our 
method as follows: first, use DMRG as done in [T3] to find the ground state on a long chain. Then, suppose we are 
interested in local observable on a region of length Iq. To find how these evolve a time t after a quench, locate some 
region of length Iq + 2vt in the chain. 

We then "observe" the state of the system outside this region to statistically sample different pure states within 
the region, and then evolve the pure states within the region. This is done as follows. The DMRG ground state can 
be written in the form 



where are a set of orthonormal states on the chain to the left of the region, |^^) are a set of orthonormal states 
to the right of the region, '^"^ are a set of states (normalized to unity but not necessarily orthogonal) on the given 
region of length Iq + 2vt, and A"^ are a set of amplitudes. We choose an a and a f3 according to the probability 



and then evolve the state 'i'^ using our present algorithm. This corresponds to observing the matrix product state 

in the basis l^*^). Repeating this statistical sampling many times we obtain the desired quantities on the local 
region. Note that on each iteration we statistically sample a,P as well as doing the sampling above. The statistical 
sampling of the state outside the region may be justified using Lieb-Robinson bounds as before. Further, when a 
matrix product state is written in the canonical form, the bond variables naturally have the orthonormality property 
used above to do the statistical sampling, giving a state in the form (|29p . 

In some cases, if the entanglement entropy grows linearly in time but with a sufhciently small prefactor, it may be 
worth using the light-cone ideas above, but doing the initial evolution for a time ti using matrix product methods 
instead of quantum circuit methods, as follows. Suppose the matrix product methods require an effort ^^2-^*, for some 
number /, to simulate for a time t. Then, to compute an observable at time tf, we simulate a subchain of length 
2vswtf for a time ti using matrix product methods. We then statistically sample states on a smaller subchain of 
length 2vsw{tf — ti) and perform the simulation of that subchain for time tf — ti exactly. The total effort is then 




(29) 



P= \A°'^\'^, 



(30) 



(31) 
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FIG. 8: Time dependence of (S'^) for the central spin as a function of time for A = 2. Exact curve is A'^ = 20 (black). Light-cone 
quantum circuit curves are I = 18 (red) and I = 20 (green). 



Choosing ti = tf/{l + /) to minimize the computational cost, we find that the cost scales as 

i^O(2^'*0, (32) 

with 

E. The Corner Transfer Quantum Circuit Method 

In this subsection we introduce the corner transfer quantum circuit method. It is primarily of theoretical, rather 
than practical interest. For calculation of local quantities (such as the expectation value of the spin on a single site) 
the light-cone method above is less work. However, the corner transfer method does give an approximation to the 
full wavefunction, and may be less work than variational matrix product methods in cases where the entanglement 
entropy grows rapidly. 

We now define the quantum circuit that approximates exp(—iHt). We define a length scale I' << VLut: the error 
in our quantum circuit approximation to exp(— ii/t/2) will be exponentially small in I', while the maximal velocity 
of information propagation for our quantum circuit will be 

Vmax = VLR + 0{l'/vLRt). (34) 

We construct the quantum circuit in two steps, first presenting a quantum circuit that approximates exjp{—iHt/2). 
We introduce operators Uk that describe the time evolution under a time dependent Hamiltonian: each Uk contains at 
a time t only the interaction terms which are contained within one of the triangles with a flattened top and jagged sides 
surrounded by a dashed line shown in in Fig. [H^a). That is, we break the time t/2 into uq = \vLii{t/2)~\ subintervals 
of time at most I/vlb.. We place the center of the triangles on sites kl, for k integer, where 



1^21' + VLRt. 



(35) 
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FIG. 9: (a) The dashed hnes show the region of space-time used in defining operators Uk, with space on the horizontal axis 
and time on the vertical axis. Only three such regions are shown, but the pattern repeats over the entire system. We also show 
this as a quantum circuit writing Ui = J7i,0j C^i,2 where each operator f/i,„ computes the exponential of the Hamiltonian 
in a given time slice. Construction is shown for I' = 2, n = 3. (b) Action of operators 14 as discussed in text, (c) Iterating 
many rounds of the corner transfer quantum circuit. On the bottom row we label U and V; on the row above, a U sits above 
each V and a V above each U. 



Then we define Uk,n for < n < no — 1 by 



i<kl+l' /2+n 

Uk,n = e^p[-i ^ hi{t/2no)], 

i=kl—l'/2—n 



(36) 



and define 



Uk = Uk,oUk,lUk,2--Uk,no-l- 



(37) 
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We then define the operators Vk as follows. We define 

kl 



C/,^„=expH E h,{t/2no)], (38) 

i—kl — I' /2 — 71 
kl+l' /2+n 



i—kl 

rLM 



That is, Uj^ represent evolution under the left or right half of the triangles in in Fig. ((HJa). Then, we define 

(fe+i)i+;/2-i ^ ^ 

Ffe=exp(-z J2 M/2)(L/,^^of/M---C^fc'!no-i)'(^^+i^o^'^+i.i---^^'^+i."o-i) ■ (39) 
i;=(fc+l/2)i-i/2+l 

That is, each Vk "undoes" the evolution in the shaded triangles as shown in Fig. Eljb), and then performs the full 
evolution in the rectangle bordered by the solid line. 
We now approximate 

exp(-^i^^/2) « n n ^'^^ (40) 

k k 

Using Lieb- Robinson bounds, one can show that the error in this approximation is 

\\exp{-iHt/2) ~l[Vkl[Uk\\ < 0{tY,\Mexp{-0{l')). (41) 

k k i 

In the second step of construction the corner transfer quantum circuit, we define another approximation to 
exp{-iHt/2). We set 

i<{k+l/2)l+l' /2+n 

Uk.n=exp[-i ^ h,{t/2no)], (42) 

i={k+l/2)l-l' /2-n 

and define 

Uk = UkfiUk,lUk,2-Uk,no-l- (43) 

Here, the centers of the upward facing triangles in Fig.[ni[a) are shifted by 1/2, compared to ([57)1 . We then define Vk 
in analogy to Eq. ([55)1 with the centers again shifted hy 1/2 and approximate 

exp(-^i/^/2) « n n ^'^^ (44) 

k k 

In Fig.[9jc) we show multiple rounds of the quantum circuit. Except for the row of triangles on the top and bottom 
boundaries, the quantum circuit looks like diamond-shaped patches in space time. They are rotated 45-degrees from 
the patches in fl\, justifying the name "corner transfer" ; the 45-degree rotation is the key improvement over [Tj leading 
to the bound on information propagation (|34l) . 



III. QUANTUM BELIEF PROPAGATION 

In this section we apply quantum belief propagation to disordered systems. We begin with a brief review of quan- 
tum belief propagation, focusing on the computational effort required, and then discuss modifications for disordered 
systems. We applied this procedure to two different disordered systems: one a chain with no frustration where Monte 
Carlo is available for comparison] 29], and the other a frustrated spin system with disorder [30|. 

Quantum belief propagation |23 is a method for constructing a matrix product density operator for a thermal 
state of a quantum system on a line or other loopless lattice. The algorithm depends on a parameter Iq and the 
approximation to the thermal state has the form in one dimension of: 

P=(o^m-i„+i-OIoI)po(o^02...0n-Io+X (45) 
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where the operators Oi act on sites i, i + Iq — 1 and the density matrix po has the form po = pl"'^° (g) 1 ig) ... (g) 1, 
with Pq "''° being a density matrix on sites l.../o- The implementation of the algorithm depends on tracing out sites, 
a process analogous to the observations discussed above. Suppose we wish to compute the partition function. If all 
of the Oi are the same, as they would be in a system without disorder, then we can consider the completely positive 
map 

P^tn{olpO,)®t), (46) 

where tri{...) denotes a trace over the first site, and p is a density operator on an interval of length Iq. We start with 
this map at p = po, and iterate it until we reach the end of the chain. We then compute the trace of the final p and 
this is the of the partition function. A similar procedure can be applied to compute expectation values. 

The cost of the procedure scales exponentially in Iq, as it requires diagonalizing matrices of dimension 2'" 26]. The 
procedure is effective down to inverse temperatures 

P - h/vLR. (47) 

The physical intuition is that the algorithm keeps quantum effects only up to a length scaling Iq. However, as we will 
see, the algorithm is capable of keeping track of classical correlations on much longer length scales; this should be no 
surprise since for classical Hamiltonians which are the sum of commuting operators, such as the Ising Hamiltonian 
H = SfSfj^i, quantum belief propagation exactly reproduces classical transfer matrix methods which can exhibit 
correlation lengths exponentially large in /3. 

For a disordered system, the operators Oi may vary as a function of i, but they depend only on the bonds on sites 
+ ^0 ~ 1- Below, we consider a system in which the bonds can assume only discrete values; in the first system 
there are 2 • 2^'°"^'/^ different possible choices for the set of bonds on sites + Iq — 1 while in the second there are 
2*0-1 pQgsibig choices. We then use the following algorithm: first, we pre-compute the operator O for each possible 
choice of bonds. We then randomly generate a configuration of bonds, and iterate the map above, choosing the 
appropriate O at each step. This requires an effort which scales linearly in system size. 

Interestingly, the algorithm seems to work better at low temperatures on disordered systems than on ordered 
systems. The reason is probably the following: when deriving the algorithm in [2^ . we used Lieb- Robinson bounds 
with velocity vlr. However, just as in the non-equilibrium case above where the actual velocity Vsw is less than VLR^ 
in these thermal systems the actual velocity may again be less than v^b.- For systems with disorder, localization 
effects may further reduce the velocity, and even change the ballistic spreading of the wavepacket to a slower growth. 
Interestingl y, t his phenomenon is known by different names in condensed matter, where it is called many- body 



localization 27|, and quantum information theory, where it is called a Lieb- Robinson bound for a disordered system [28 



A. Results on Disordered Systems 



The first disordered system we considered has the Hamiltonian 

N/2 N/2 

H = J S2i-i ■ S2i + JiS2i ■ S2i+i, (48) 

1=1 2=1 

where each spin is spin- 1/2, J > and Ji = Jp < Q with probability p and Ji — Ja > ^ with probability 1 — p. 
This model has both ferromagnetic and anti-ferromagnetic couplings, and was proposed to model the compound 
{CH^)2CHNH^Cu{ClxBri^x)?,i where the probability p = a;^ [H, Hj, [H, [H, [s^l ■ The case considered is J = 1 and 
\Ji\ = 2. We performed simulations on this chain with Iq = 5, 7, 9 for p = 0.2, 0.4, 0.6, 0.8 and for Iq = 5, 7, 9, 11 for 
p = 0, 1. For the random cases {p — 0.2, 0.4, 0.6, 0.8) we considered chains of 100000 sites and computed the uniform 
susceptibility by the change in partition function in response to a weak applied field. In this way, the susceptibility 
was self-averaging. For the pure cases, we considered shorter chains, and computed the applied susceptibility by 
measuring partition functions as in The results are shown in Figs. llOiTT] For small Iq, qualitatively wrong results 
are seen at low temperature, with the susceptibility diverging at low enough temperature. However, as is increased, 
the accuracy extends to lower temperature. The difference between the curves for — 7 and lg=9 is small for T 
greater than roughly 1/7. In this region also, we agree well with Monte Carlo data. 

Data was taken over several different (3, with a step of 0.25 in /3. As in [g^l, we set the perturbation A (following 
notation of [l^) to equal (l/2)(/i;„_2. /n-i + /i;o-i,io) rather than A ~ /ijo-i,io- 

The second model we considered is [30] a model with frustration and disorder, where quantum Monte Carlo results 
are not available. The study in [30J was motivated by experimental studies on CuGeO^^^ and CuQGeeOis — xH20 3'^, 
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FIG. 10: Uniform susceptibility as a function of temperature for different p,lo. Top: lo = 5 (black), lo = 7 (red), lo ~ 9 (green). 
Middle: lo — 5 (black), lo — 7 (red), lo — 9 (green), Zq = 11 (blue). Bottom: lo = 5 (black), lo — 7 (red), lo — 9 (green), lo — 11 
(blue). 



where second neighbor interactions may be important. We were interested in a case where second neighbor interactions 
would be very important, so we considered the Hamiltonian 

JV-l JV-2 

H=Y, J^S, ■ S,+i + K^S, ■ S,+2, (49) 
1=1 1=1 

and in the pure case we considered Ji — l,Ki = 1/2. This is a Majumdar-Ghosh chain with a dimerized ground 
state. In the disordered case we chose Ji — 0.9 or J; = 1.1, with probabihty 1/2 of either choice (similar results were 
found for choosing J = 0.75 or J = 1.25). We choose 

= (l/2)J,J,+i, (50) 

which correlates the second neighbor interaction with the nearest neighbor interaction as described in [2^ . 

We studied Iq — 5, 7,9 with chains of length 19999. It is necessary to take such long chains in the pure case to 
avoid boundary condition effects because at low temperatures in the pure system there is an exponentially increasing 
correlation length for dimer-dimer correlations. In Fig. [T2t we first show the results of the specific heat as a function 
of /3 for the pure case, computed from the second derivative of the partition function (probably a very slightly more 
accurate method is to take the first derivative of the energy as in [22|). A strong difference is seen between Iq — 5 and 
Iq — 7 above /3 ^ 3.25, but the Iq — 7 and Iq — % curves are almost identical. This indicates that by going to Zq = 9 
we have succeeded in converging the specific heat in Iq for (3 < 10. 

The uniform susceptibility shows a similar effect. The pure curves show a large difference between Iq — 5 and 
Iq = 7, 9, but only slight difference can be seen between Zq = 7, 9 and only above (3 = 8. Again, the results seem to be 
converged in Iq by going to Zq = 9 in the range of temperatures we consider. The disordered curves show again that 
Iq = 5 is too small, but for Iq = 7,9 little difference is seen (except for some small random fluctuations) up to /3 = 10. 
A very slightly higher susceptibility is seen in the random case compared to the pure case. Finally, we consider the 
dimer susceptibility, defined by 

N/2-1 ^ ^ 2 

Xdimer — S2i ■ S2i+1 — 821+1 ■ 821+2^ )■ (51) 

i=l 
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FIG. 11: Uniform susceptibility as a function of temperature for p — 0.2,0.4,0.6 from top to bottom and for Iq = 5 (black), 
lo ^7 (red), Zo = 9 (green). 



This shows a large difference between the pure and disordered cases. The pure cases again show agreement between 
^0 = 7, 9 and show an an increase that gives Xdimer/ P growing exponentially in (3. The disordered cases show Xdimer/P 
saturating as a function of (3. 

The saturation of the dimer-dimer correlation function in the disordered case is no surprise. The disorder locally 
breaks the Z2 symmetry between different ordering patterns of the dimers, and is relevant for this one dimensional 
system. The fact that the uniform susceptibility shows only very slight difference between pure and disordered systems 
is more surprising; asymptotically, the uniform susceptibility should decay exponentially in (3 in the pure system and 
should increase as /3/log^(/3) in the random singlet phase [30|. 



IV. DISCUSSION 



The main result in this paper is the light-cone quantum circuit method. We have tested this method numerically 
on a free system, with A = 0, and on interacting systems with A ^ 0. We have found decaying oscillations in the 
expectation value of the spin. In future, this technique will be useful for studying non-integrable systems, to see if 
they relax to a thermal state [sst. 

We can study two-dimensional systems by considering them as wide one-dimensional systems; this allows us to 
double the number of spins, but only leads to a factor of -\/2 increase in the time compared to direct simulation. 
Other similar quantum circuit methods may be more effective in two dimensions. 

The results using the light-cone quantum circuit method are indeed comparable to those one would find by exactly 
solving a system of twice the size, at least for A = where we can compute exactly. Thus, we I = 18, we find results 
comparable to a system of TV = 35 or TV = 37 sites. The computational cost to exactly evolve a given system is roughly 
comparable to that required to do exact diagonalization using Lanczos methods on that system: Lanczos methods and 
exact evolution both require sparse matrix-vector multiplication, but the number of multiplications needed to reach 
convergence for the time evolution may be larger than that needed to reach convergence for ground state properties. 
Thus, we expect that sizes around 35 sites, especially given the low symmetry of the present system, are around the 
upper limit for exact methods now, while we carried the light-cone quantum circuit method up to Z = 22. Further, the 
asymptotic analysis of time requirements applies also to memory requirements: the memory requirements of an exact 
solution on N sites scale as N2^ , while while the memory requirement of the light-cone quantum circuit method scale 
as 12\ so, regardless of what N can be obtained using an exact solution, it should be possible to obtain the same /, up 
to a difference of a couple sites, in the light-cone quantum circuit method. The main additional cost in the light-cone 
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FIG. 12: Top: specific heat for the pure system. Middle: uniform susceptibility for the pure system {Iq = 5, 7, 9 are black, red, 
green respectively) and the disordered system {Iq — 5,7,9 are blue, yellow, brown respectively). Bottom: dimer susceptibility 
for the pure system {lo = 5, 7, 9 are black, red, green respectively) and the disordered system {Iq = 5, 7, 9 are blue, yellow, 
brown respectively). 



quantum circuit method is the need to run many times to obtain statistical samples, but this is a problem which can 
be parallelized. 
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